NumEx2

In this numerical exercise, you will use the galaxy and particle/cell information from the cosmological simulation IllustrisTNG. In particular, wee will use IllustrisTNG100-1 and IllustrisTNG50-1 (hereafter TNG100, and TNG50). The numerical exercise is divided into two activities. In the first part, you will evaluate and test different scaling relations from the galaxy catalogs with observational results. In the second part, we will play with the information from the particles to investigate the inner structure of the different components.

First steps: Log in to the IllustrisTNG database (https://www.tng-project.org/, Springel et al., 2018) and enter the JupyterLab application. Upload the scripts from our sha- ring folder into your JupyterLab. The scripts will be the basis for our analysis. Through the JupyterLab you can obtain information on dark-matter haloes, galaxies, and particle infor- mation at different redshifts (snapshots). We will start using snapshot 99 which represents z = 0.

Activity 1

1. Compute the present-day (i.e., z = 0) Galaxy Stellar Mass Function for TNG50-1 and TNG100-1. Compare your results with the observations from Figure 1 (you can use: https://apps.automeris.io/wpd/ for extracting the data).

#####################################
import warnings
warnings.filterwarnings("ignore")

import logging
logging.getLogger().setLevel(logging.CRITICAL)

PRINT=True

import h5py
import numpy as np
from numpy import linalg as LA 
import matplotlib.pylab as plt
from matplotlib.colors import LogNorm
import illustris_python as il
import matplotlib.gridspec as gridspec
import matplotlib as mpl
from scipy.stats import median_abs_deviation as MAD # Median ABbsolute Deviation


%matplotlib inline

plt.style.use('dark_background')
plt.style.use('seaborn-bright')

COLOR = 'lightgray'

mpl.rcParams['text.color']      = COLOR  # labels, titles, etc colors
mpl.rcParams['axes.labelcolor'] = COLOR  # Numbers in axis colors
mpl.rcParams['xtick.color']     = COLOR  # ticks colors
mpl.rcParams['ytick.color']     = COLOR  # ticks colors
mpl.rc(      'axes',edgecolor   = COLOR) # axis Colors

mpl.rc('font', size=15) # Font Size

mpl.rcParams['figure.figsize'] = [10, 10] # Size of figures

Data={}
LLL =[50,100,300]

VVV=[51.7,110.7,302.6]
h  = 0.6774


POSITION= ['Fisrt','Second','Third']

# LLL =[50]
for i in LLL:
    if PRINT:
        print(f'Loading data of TNG{i}-1 ...   ')
    basePath = '/home/tnguser/sims.TNG/TNG'+str(i)+'-1/output/'
    snapNum  = 99  # z=0
    fields   = ['SubhaloFlag', 'SubhaloHalfmassRadType', 'SubhaloMassType', 'SubhaloPos', 'SubhaloVel',\
                'SubhaloSFRinRad', 'SubhaloSFR','SubhaloStarMetallicity']
    
    
    Subhalo  = il.groupcat.loadSubhalos(basePath=basePath,snapNum=snapNum,fields=fields)
    
    # SubhaloFlag = 0 means no cosmological origin
    # stellar mass range log(M/h) = 9.0 to above
    
    mask1= (Subhalo['SubhaloFlag']!=0)
    
    idx = []
    R200= []
    M200= []
    DM  = Subhalo['SubhaloMassType'][:,1].copy()[mask1]
    for PPP in range(3):
        IDX=DM.argmax()
        idx.append(IDX)
        if PRINT:
            print(f'In TNG{i}-1 the {POSITION[PPP]}\t most massive Halo ({IDX:^10}) has:\t {np.log10(DM[IDX]/h)+10:.3f}  [Log(M_sun)] ')
        DM[IDX]= 0 
        
    del DM

    
    mask = np.where((Subhalo['SubhaloMassType'][:,4]>=0.003) & mask1)
    
    globals()['mass_st_'+str(i)] = Subhalo['SubhaloMassType'][:,4][mask]/h # to save the mass in units of 10^10 solar masses 
    
    Data[str(i)]= dict( sh_ids   = mask[0], # The subhalo "ID" is the same as the "index"
                        cm       = Subhalo['SubhaloPos'][mask],
                        cv       = Subhalo['SubhaloVel'][mask],
                        rh_type  = Subhalo['SubhaloHalfmassRadType'][:,4][mask],
                        sfr_2rh  = Subhalo['SubhaloSFRinRad'][mask],     # in 2 * R_half
                        sfr_tot  = Subhalo['SubhaloSFR'][mask],          # in 2 * R_half
                        mh       = Subhalo['SubhaloStarMetallicity'][mask],
                        MassiveH = np.array(idx))

    
    # clear space
    del Subhalo
    if PRINT:
        print('Done')
        print()

Loading data of TNG50-1 ...   
In TNG50-1 the Fisrt	 most massive Halo (    0     ) has:	 14.225  [Log(M_sun)] 
In TNG50-1 the Second	 most massive Halo (  63784   ) has:	 13.897  [Log(M_sun)] 
In TNG50-1 the Third	 most massive Halo (  96666   ) has:	 13.757  [Log(M_sun)] 
Done

Loading data of TNG100-1 ...   
In TNG100-1 the Fisrt	 most massive Halo (  17066   ) has:	 14.545  [Log(M_sun)] 
In TNG100-1 the Second	 most massive Halo (    0     ) has:	 14.537  [Log(M_sun)] 
In TNG100-1 the Third	 most massive Halo (  31116   ) has:	 14.464  [Log(M_sun)] 
Done

Loading data of TNG300-1 ...   
In TNG300-1 the Fisrt	 most massive Halo (    0     ) has:	 15.205  [Log(M_sun)] 
In TNG300-1 the Second	 most massive Halo (  11704   ) has:	 15.047  [Log(M_sun)] 
In TNG300-1 the Third	 most massive Halo (  17827   ) has:	 14.974  [Log(M_sun)] 
Done

def mid(x):
    return (x[1:]+x[:-1])/2

# From the graph
D1,D2 = np.genfromtxt('full_line.csv',delimiter=',',unpack=True)

l,u  = np.log10(min(D1)),np.log10(max(D1)) # Range of masses desired 

G1,G2 = np.genfromtxt('dots.dat',delimiter=',',unpack=True)

# From the NUMEX1
F1,F2,F3,F4,F5,F6 = np.genfromtxt('MyHMF.dat',delimiter='\t',unpack=True)


###################################################################
# plt.plot(D1,D2,label='Ilustris_TNG_paper',marker='.',color='k')

# plt.legend()
# plt.xscale('log')
# plt.yscale('log')

# plt.show(),plt.close()

# print(l,u)

fig= plt.figure(figsize=(18,6))
gs = gridspec.GridSpec( 1,5,                             # Number of axis y,x
                        height_ratios = [1],             # relatives ratios of heigh
                        width_ratios  = [1,0.1,1,0.1,1], # relatives ratio of with
                        left  = 0.1,      # Space to the edge of left   from the nearest axis
                        right = 0.95,     # Space to the edge of right  from the nearest axis
                        bottom= 0.15,     # Space to the edge of bottom from the nearest axis
                        top   = 0.97,     # Space to the edge of top    from the nearest axis
                        wspace= 0.1,      # Space horizontal between each of the axis
                        hspace= 0.0)      # Space vertical   between each of the axis

ax1  = fig.add_subplot(gs[0,0])
ax2  = fig.add_subplot(gs[0,2],sharey=ax1,sharex=ax1)
ax3  = fig.add_subplot(gs[0,4],sharey=ax1,sharex=ax1)

AX=[ax1,ax2,ax3]

# AX[0].ticklabel_format(axis='y', scilimits=[-1, 1])
AX[0].set_ylim(1e-5,0.1)

# Ploting the masses
for j,i in enumerate(LLL):
    VV  = (VVV[j])**3
    XX  = np.log10(globals()['mass_st_'+str(i)])+10##np.append(log_mass_st_q,log_mass_st_sf)
    XX  = XX[(XX>l)&(XX<u)]
    y,x = np.histogram(XX,bins=50)
    w   = x[1]-x[0]
    
    AX[j].plot(10**mid(x),y/VV/w,label=f'TNG {i}-1',marker='s') 
    
    AX[j].plot(D1,D2,label='Ilustris_TNG_paper',marker='.',color='gray')
    AX[j].scatter(G1,G2,label='GAMA-II (Wright et al. 2017)',marker='o',edgecolor='r',facecolor='none')
    
    
    AX[j].legend(loc='upper left', bbox_to_anchor=(0.27, 1.12),fontsize=12)
    
    AX[j].grid(zorder=0,ls=':',color='gray')
    AX[j].set_yscale('log')
    AX[j].set_xscale('log')

fig.supylabel(r'$N_{Halo(M_Stars)}\cdot\ dv^{-1}\cdot\ d(Log_{10}(M))^{-1}$',fontsize=18)
fig.supxlabel('Log (Stellar Mass [M$_{\odot}$])',fontsize=18)
print()
plt.show()
print('\n\n\n')

2. Compare your results with the halo mass function obtained in Numerical Exercise 1. What differences do you see in the GSMF and halo mass function distributions? Interpret the results of the GSMF shape in terms of the physical processes that regulate the stellar mass of the galaxies.

fig= plt.figure(figsize=(17,5))
gs = gridspec.GridSpec( 1,5,                             # Number of axis y,x
                        height_ratios = [1],             # relatives ratios of heigh
                        width_ratios  = [1,0.1,1,0.1,1], # relatives ratio of with
                        left  = 0.08,      # Space to the edge of left   from the nearest axis
                        right = 0.95,     # Space to the edge of right  from the nearest axis
                        bottom= 0.18,     # Space to the edge of bottom from the nearest axis
                        top   = 0.97,     # Space to the edge of top    from the nearest axis
                        wspace= 0.1,      # Space horizontal between each of the axis
                        hspace= 0.0)      # Space vertical   between each of the axis

ax1  = fig.add_subplot(gs[0,0])
ax2  = fig.add_subplot(gs[0,2],sharey=ax1,sharex=ax1)
ax3  = fig.add_subplot(gs[0,4],sharey=ax1,sharex=ax1)

AX=[ax1,ax2,ax3]

# AX[0].ticklabel_format(axis='y', scilimits=[-1, 1])
AX[0].set_ylim(1e-5,1)

# Ploting the masses
for j,i in enumerate(LLL):
    VV  = (i)**3
    XX  = np.log10(globals()['mass_st_'+str(i)])+10##np.append(log_mass_st_q,log_mass_st_sf)
    XX  = XX[(XX>l)&(XX<u)]
    y,x = np.histogram(XX,bins=50)
    w   = x[1]-x[0]
    
    AX[j].plot(10**mid(x),y/VV/w,label=f'TNG {i}-1 (Stars)',marker='s',zorder=10) 
    
    AX[j].plot(10**F1,F2,label='SMP (Halo)',marker='.')
    AX[j].plot(10**F3,F4,label='TNG 300 (Halo)',marker='.')
    AX[j].plot(10**F5,F6,label='EAGLE (Halo)',marker='.')
    
    
    AX[j].legend(loc='upper right',fontsize=12)#, bbox_to_anchor=(0.52, 1.12))
    
    AX[j].grid(zorder=0,ls=':',color='gray')
    AX[j].set_yscale('log')
    AX[j].set_xscale('log')

fig.supylabel(r'$N_{Halo/Stars}\cdot\ dv^{-1}\cdot\ d(Log_{10}(M))^{-1}$',fontsize=18)
fig.supxlabel('Log (Mass [M$_{\odot}$])',fontsize=18)
print()

first we can see an obvious offset of the mass, and the other characteristic that we can see is a break around $3\times10^{10} M_{\odot}$ where the distribution for lower masses are dominated by the Super nova feedback while for the larger masses are more affected by AGN feedback

print('\n\n\n')

3. We now focus on different scaling relations of the galaxies. Map the stellar mass vs. star formation rate of the galaxies in TNG50 and TNG100. Identify the star-forming and quenched galaxies using a threshold in $sSFR = 10^{−11} yr^{−1}$. You can use the script available.

np.random.seed(1209)

Save_Galaxies={}

sub_data={}
for i in Data.keys():
    # starformin and quenched subhalos [https://arxiv.org/pdf/1812.07584.pdf]
    mass_st      = globals()['mass_st_'+str(i)]
    MH           = Data[i]['mh']
    sfr_tot      = Data[i]['sfr_tot']
    log_ssfr_tot = np.log10(sfr_tot / mass_st / 10**10)
    mask_sf      = (log_ssfr_tot > -11)
    mask_q       = (log_ssfr_tot < -11)
    sub_data[i]=dict(log_mass_st_sf = np.log10(mass_st[mask_sf]) + 10,
                     log_mass_st_q  = np.log10(mass_st[mask_q])  + 10,
                     log_sfr_tot_sf = np.log10(sfr_tot[mask_sf]),
                     log_sfr_tot_q  = np.log10(sfr_tot[mask_q]),
                     MH_sf          = MH[mask_sf],
                     MH_q           = MH[mask_q])
    
    # This is for activity 2.1
    ids       = Data[i]['sh_ids']
    mask_mass = (-1<np.log10(mass_st))&(np.log10(mass_st)<1) 
    
    
    SF        = np.random.choice(ids[(mask_sf & mask_mass)],10)
    Q         = np.random.choice(ids[(mask_q  & mask_mass)],10)
    
    Save_Galaxies[str(i)]={'SF':SF,'Q':Q}
    
    if PRINT:
        print(f'TNG {i}\nStar Forming: \t{SF}\nQuenched:  \t{Q}')
        print()

TNG 50
Star Forming: 	[632978 720707 386271 419620 627888 708395 167435 733149 300910 562743]
Quenched:  	[504559     46 282781 470348  63937  96776 117259 289398  63875 143904]

TNG 100
Star Forming: 	[577306 533722 419303 630829 595123 606391 580474 570763 332070 337846]
Quenched:  	[134867  76208 279723 496713 182005  52859 114435 458028 197149 366793]

TNG 300
Star Forming: 	[2116125 2168130 2511834 2256571  895130 1167096 1940368 1741389 2302679
 1249631]
Quenched:  	[ 902307   92281  834988 1232729  786460  904851  264577  777808  299342
 1425898]

fig= plt.figure(figsize=(17,5))
gs = gridspec.GridSpec( 1,5,                             # Number of axis y,x
                        height_ratios = [1],             # relatives ratios of heigh
                        width_ratios  = [1,0.1,1,0.1,1], # relatives ratio of with
                        left  = 0.08,      # Space to the edge of left   from the nearest axis
                        right = 0.95,     # Space to the edge of right  from the nearest axis
                        bottom= 0.18,     # Space to the edge of bottom from the nearest axis
                        top   = 0.97,     # Space to the edge of top    from the nearest axis
                        wspace= 0.1,      # Space horizontal between each of the axis
                        hspace= 0.0)      # Space vertical   between each of the axis

ax1  = fig.add_subplot(gs[0,0])
ax2  = fig.add_subplot(gs[0,2],sharey=ax1,sharex=ax1)
ax3  = fig.add_subplot(gs[0,4],sharey=ax1,sharex=ax1)

AX=[ax1,ax2,ax3]

# AX[0].ticklabel_format(axis='y', scilimits=[-1, 1])
#AX[0].set_ylim(1e-5,1)

Var= ["log_mass_st_sf" ,"log_mass_st_q" , "log_sfr_tot_sf" ,"log_sfr_tot_q","MH_sf","MH_q"]

# Ploting the masses
for j,i in enumerate(LLL):
    
    DDD = sub_data[str(i)]
    
    AX[j].scatter(DDD[Var[0]], DDD[Var[2]], s=1, color='c', label='Star-Forming')
    AX[j].scatter(DDD[Var[1]], DDD[Var[3]], s=1, color='m', label='Quenched')
    
    AX[j].legend(loc='upper left',title=f'TNG{i}-1',fontsize=12, bbox_to_anchor=(.0, 1.12),framealpha=1)
    
    AX[j].grid(zorder=0,ls=':',color='gray')
    # AX[j].set_yscale('log')
    # AX[j].set_xscale('log')
          

fig.supylabel('Log (SFR [M$_{\odot}$/yr])',fontsize=18)
fig.supxlabel('Log (Stellar Mass [M$_{\odot}$])',fontsize=18)
print()
plt.show()
print('\n\n\n')

4. For each sub-set, plot: i) the half-mass stellar radius vs. stellar mass and compute the median of each subset and simulation; ii) the mass-metallicity relation for star-forming vs. quenched galaxies. Interpret the results in each case.

fig= plt.figure(figsize=(18,9))
gs = gridspec.GridSpec( 3,3,                            # Number of axis y,x
                        height_ratios = [0.5,0.04,0.5],  # relatives ratios of heigh
                        width_ratios  = [0.33,0.33,0.33],               # relatives ratio of with
                        left  = 0.1,      # Space to the edge of left   from the nearest axis
                        right = 0.95,     # Space to the edge of right  from the nearest axis
                        bottom= 0.1,     # Space to the edge of bottom from the nearest axis
                        top   = 0.97,     # Space to the edge of top    from the nearest axis
                        wspace= 0.2,      # Space horizontal between each of the axis
                        hspace= 0.0)      # Space vertical   between each of the axis

ax1  = fig.add_subplot(gs[0,0])
ax2  = fig.add_subplot(gs[2,0])
ax3  = fig.add_subplot(gs[0,1],sharex=ax1,sharey=ax1)
ax4  = fig.add_subplot(gs[2,1],sharex=ax2,sharey=ax2)
ax5  = fig.add_subplot(gs[0,2],sharex=ax1,sharey=ax1)
ax6  = fig.add_subplot(gs[2,2],sharex=ax2,sharey=ax2)

# (here we can define if we want to share axes and with which one)

AX=[ax1,ax2,ax3,ax4,ax5,ax6]

NNN= 9

#--
# Ploting
Var2=["sh_ids","cm","cv","rh_type","sfr_2rh","sfr_tot","mh"]

if PRINT:
    print('#--------------------------------------------------------------')
    print('#SET\tMedian Mass\t Median HMR\t MAD mass\t MAD HMR')
    print('#   \t Log(mass) \t Log(HMR)  \tLog(mass)\t Log(HMR)')


for j,i in enumerate(LLL):
    j=j*2
    mass_st = np.log10(globals()['mass_st_'+str(i)])+10
    HMR     = np.log10(Data[str(i)][Var2[3]])
    
    u,l= max(mass_st),min(mass_st)
 
    AX[j].scatter(mass_st, HMR, s=1, color='b')#, label='HMR')
    # AX[j].legend(loc='upper left',fontsize=12, bbox_to_anchor=(.0, .98),framealpha=1)
    Lsp  = np.linspace(l,u+0.1,NNN+1)
    for O in range(NNN):
        
        mask = (Lsp[O]<=mass_st)&(mass_st<Lsp[O+1])
        
        x,y   = np.median(mass_st[mask]),np.median(HMR[mask])
        XE,YE = MAD(mass_st[mask]),MAD(HMR[mask])
        
        if PRINT:
            print(f'TNG{i:^4}: {x:.4f}   \t {y:.4f}   \t {XE:.4f}  \t {YE:.4f}')
    
        AX[j].errorbar(x,y, yerr=YE, xerr=XE,color='r')
    
    DDD = sub_data[str(i)]
    
    AX[j+1].scatter(DDD[Var[0]], np.log10(DDD[Var[4]]), s=1, color='c', label='Star Forming',alpha=0.5)
    AX[j+1].scatter(DDD[Var[1]], np.log10(DDD[Var[5]]), s=1, color='m', label='Quenched',alpha=0.5)
    
    AX[j].set_title(f'TNG{i}-1')
    AX[j+1].legend(loc='lower right',fontsize=12, bbox_to_anchor=(1, 0),framealpha=1)
    
    AX[j].grid(zorder=0,ls=':',color='gray')
    AX[j+1].grid(zorder=0,ls=':',color='gray')
    # AX[j].set_yscale('log')
    # AX[j].set_xscale('log')
    
    if PRINT: 
        print()
    
    
AX[1].set_ylim(-3.2,-1)    

    
# fig.supylabel('Log (SFR [M$_{\odot}$/yr])',fontsize=18)

fig.supxlabel('Log (Stellar Mass [M$_{\odot}$])',fontsize=18)

# Press and Schechter


AX[0].set_ylabel("half-mass stellar radius Log(cKpc/h)",labelpad=12)
AX[1].set_ylabel("Metallicity  Log(Z)",labelpad=12)


# Erasing extra y labels
plt.setp(ax1.get_xticklabels(), visible=False)
plt.setp(ax3.get_xticklabels(), visible=False)
plt.setp(ax5.get_xticklabels(), visible=False)
print()
plt.show()
print('\n\n\n')

#--------------------------------------------------------------
#SET	Median Mass	 Median HMR	 MAD mass	 MAD HMR
#   	 Log(mass) 	 Log(HMR)  	Log(mass)	 Log(HMR)
TNG 50 : 7.9083   	 -0.1485   	 0.1402  	 0.1143
TNG 50 : 8.4990   	 -0.0979   	 0.1400  	 0.1451
TNG 50 : 9.0477   	 0.0990   	 0.1330  	 0.1645
TNG 50 : 9.6165   	 0.3084   	 0.1446  	 0.1763
TNG 50 : 10.2156   	 0.4381   	 0.1384  	 0.2115
TNG 50 : 10.7552   	 0.5779   	 0.1398  	 0.1613
TNG 50 : 11.2716   	 0.7661   	 0.1269  	 0.1085
TNG 50 : 11.8787   	 0.9771   	 0.1021  	 0.0907
TNG 50 : 12.7353   	 1.4666   	 0.0000  	 0.0000

TNG100 : 7.9000   	 0.3181   	 0.1392  	 0.0949
TNG100 : 8.4642   	 0.2740   	 0.1369  	 0.0852
TNG100 : 9.0157   	 0.2729   	 0.1358  	 0.0914
TNG100 : 9.5586   	 0.3215   	 0.1336  	 0.0980
TNG100 : 10.1140   	 0.3576   	 0.1354  	 0.1250
TNG100 : 10.6488   	 0.4736   	 0.1202  	 0.1585
TNG100 : 11.1681   	 0.8514   	 0.1179  	 0.1022
TNG100 : 11.7347   	 1.1998   	 0.1123  	 0.0862
TNG100 : 12.2431   	 1.5139   	 0.1043  	 0.1265

TNG300 : 7.8887   	 0.6524   	 0.1417  	 0.1237
TNG300 : 8.5125   	 0.5498   	 0.1491  	 0.1060
TNG300 : 9.1031   	 0.4868   	 0.1455  	 0.0839
TNG300 : 9.7022   	 0.4821   	 0.1470  	 0.0857
TNG300 : 10.2998   	 0.4918   	 0.1494  	 0.0781
TNG300 : 10.7860   	 0.6862   	 0.1098  	 0.1029
TNG300 : 11.3803   	 1.0967   	 0.1124  	 0.1074
TNG300 : 11.9668   	 1.5120   	 0.0949  	 0.1052
TNG300 : 12.5110   	 1.8852   	 0.0867  	 0.1512

Activity 2

We will now use the information from the gas cells and stellar and dm particles.

1. Using the script provided, select a subset of 20 galaxies (10 star-forming and 10 quenched) with stellar mass in the range of $10^{9} − 10^{11} M_\odot$ and plot the galaxies face-on. Use the script to see the distribution of gas and stars.

if PRINT:
    print(Save_Galaxies['50'])

{'SF': array([632978, 720707, 386271, 419620, 627888, 708395, 167435, 733149,
       300910, 562743]), 'Q': array([504559,     46, 282781, 470348,  63937,  96776, 117259, 289398,
        63875, 143904])}

# Some relevant functions
def get_basic_info_of_subhalo(Sim='50', ID=None):
    DDD = Data[Sim]
    if (Sim==None or ID==None):
        print('please enter the index and ID of subhalo')
    
    mask= np.where(DDD['sh_ids']==ID)[0][0]
    
    subhalo_id = ID
    r_half     = DDD['rh_type'][mask]  # half stellar mass radius
    com        = DDD['cm'][mask]
    cov        = DDD['cv'][mask]
    return subhalo_id, r_half, com, cov

def cal_specific_ang_mom(mass, pos, vel, r_out):
    r    = LA.norm(pos, axis=1)
    mask = np.where(r <= r_out)[0]
    
    m = mass[mask]
    p = pos[mask]
    v = vel[mask]
    L = np.cross(p, v)
    L_total = np.sum(m * L.T, axis=1)
    return L_total

def vector_rotation_matrix(vec1, vec2):
    # unit vectors
    v1 = vec1 / LA.norm(vec1)
    v2 = vec2 / LA.norm(vec2)
    # dimension and identity
    dim = v1.size
    I = np.identity(dim)
    # cos of angle between v1 and v2
    c = np.dot(v1, v2)
    # cross product matrix of a vector to rotate around
    K = np.outer(v2, v1) - np.outer(v1, v2)
    # Rodrigues' formula
    return I + K + (K @ K) / (1 + c)

i=50
basePath = '/home/tnguser/sims.TNG/TNG'+str(i)+'-1/output/'

Quen_pos={}
Quen_vel={}
Quen_posG={}
Quen_velG={}

DM_or_gas_Q={}

MASS_DM = ((3.64755609660833e-05)/h)*(1e10) # https://www.tng-project.org/data/downloads/TNG50-1-Dark/

if PRINT:
    print(MASS_DM)

for i in Save_Galaxies['50']['Q']:
    fig= plt.figure(figsize=(18,9))
    gs = gridspec.GridSpec( 3,7,                             # Number of axis y,x
                            height_ratios = [1,0.23,1],             # relatives ratios of heigh
                            width_ratios  = [1,0.1,1,0.1,1,0.1,0.08], # relatives ratio of with
                            left  = 0.1,      # Space to the edge of left   from the nearest axis
                            right = 0.95,     # Space to the edge of right  from the nearest axis
                            bottom= 0.15,     # Space to the edge of bottom from the nearest axis
                            top   = 0.97,     # Space to the edge of top    from the nearest axis
                            wspace= 0.1,      # Space horizontal between each of the axis
                            hspace= 0.0)      # Space vertical   between each of the axis

    ax1  = fig.add_subplot(gs[0,0])
    ax2  = fig.add_subplot(gs[0,2])
    ax3  = fig.add_subplot(gs[0,4])
    cax1 = fig.add_subplot(gs[0,6])
    
    ax4  = fig.add_subplot(gs[2,0])
    ax5  = fig.add_subplot(gs[2,2])
    ax6  = fig.add_subplot(gs[2,4])
    cax2 = fig.add_subplot(gs[2,6])
    

    axs=[ax1,ax2,ax3,ax4,ax5,ax6]
    
    
    fig.suptitle(f'Galaxy {i} / Quenched',y=1.05)
    #################################################################################
    
    subhalo_id, r_half, com, cov = get_basic_info_of_subhalo(ID=i)
    fields = ['Coordinates', 'Velocities', 'Masses']

    st = il.snapshot.loadSubhalo(basePath,snapNum,subhalo_id,'stars',fields=fields)

    mass = st['Masses']
    pos  = st['Coordinates']
    vel  = st['Velocities']

    pos = pos - com  # shifting co-ordinate system to center of subhalo
    vel = vel - cov


    # plotting three projection of subhalo after rotation
    #fist rotate subhalo 
    L_vec = cal_specific_ang_mom(mass, pos, vel, 10*r_half)  # angular momentum vector
    rot_mat = vector_rotation_matrix(L_vec, [0,0,1])  # [0,0,1] for face-on
    new_pos = np.dot(rot_mat, pos.T)
    new_vel = np.dot(rot_mat, vel.T)
    new_pos, new_vel = new_pos.T, new_vel.T
    
    Quen_pos[str(i)]=new_pos
    Quen_vel[str(i)]=new_vel
    
        
    extent = 5  # in unit of half mass radius
    rng = [[-extent*r_half,extent*r_half],[-extent*r_half,extent*r_half]]
    nbins = int((rng[0][1] - rng[0][0])/0.12)
    
    H1,X1,Y1= np.histogram2d(new_pos[:,0], new_pos[:,1], bins=nbins, range=rng)

    H1[H1<=0]=np.nan

    axs[0].imshow(np.log10(H1.T), vmin=np.log10(3) ,cmap='gist_ncar',interpolation = 'spline36',interpolation_stage='rgba',origin='lower',extent=[-extent*r_half,extent*r_half,-extent*r_half,extent*r_half])
    axs[0].set_xlabel('x [kpc]')
    axs[0].set_ylabel('y [kpc]')
    
    
    H1,X1,Y1= np.histogram2d(new_pos[:,0], new_pos[:,2], bins=nbins, range=rng)

    H1[H1<=0]=np.nan

    axs[1].imshow(np.log10(H1.T), vmin=np.log10(3) ,cmap='gist_ncar',interpolation = 'spline36',interpolation_stage='rgba',origin='lower',extent=[-extent*r_half,extent*r_half,-extent*r_half,extent*r_half])
    axs[1].set_xlabel('x [kpc]')
    axs[1].set_ylabel('z [kpc]')
    
    H1,X1,Y1= np.histogram2d(new_pos[:,1], new_pos[:,2], bins=nbins, range=rng)

    H1[H1<=0]=np.nan

    IM=axs[2].imshow(np.log10(H1.T), vmin=np.log10(3) ,cmap='gist_ncar',interpolation = 'spline36',interpolation_stage='rgba',origin='lower',extent=[-extent*r_half,extent*r_half,-extent*r_half,extent*r_half])
    axs[2].set_xlabel('y [kpc]')
    axs[2].set_ylabel('z [kpc]')
    
    plt.colorbar(IM,cax=cax1,ticks=[np.log10(3)*1.05,np.log10(np.nanmax(H1))*0.95])
    cax1.set_ylabel('Density', rotation=270)
    cax1.set_title('Stars')
    cax1.set_yticklabels(['low','high'])
    
    ####################################################################################
    ####################################################################################
    ##############    GAS    ###########################################################
    ####################################################################################
    ####################################################################################
    
    
    try:
        st = il.snapshot.loadSubhalo(basePath,snapNum,subhalo_id,'gas',fields=fields)
        
        if st['count']<5000:
            raise Exception('To few particles showing instead DM') 
        
        mass = st['Masses']
        pos  = st['Coordinates']
        vel  = st['Velocities']
        cax2.set_title('Gas')
        
        print('Gas: Count: ',st['count'])
        
        CMAP = 'inferno'
        
        DM_or_gas_Q[str(i)]='Gas'
    except Exception as e:
        print(e)
        
        print('Gas: Count: ',st['count'])
        
        fields = ['Coordinates', 'Velocities']
        st = il.snapshot.loadSubhalo(basePath,snapNum,subhalo_id,'dm',fields=fields)
    
        pos  = st['Coordinates']
        vel  = st['Velocities']
        mass = np.repeat(MASS_DM,len(pos))
        cax2.set_title('DM')
        CMAP = 'viridis'
        
        DM_or_gas_Q[str(i)]='Dark Matter'

    pos = pos - com  # shifting co-ordinate system to center of subhalo
    vel = vel - cov


    # plotting three projection of subhalo after rotation
    #fist rotate subhalo 
    L_vec = cal_specific_ang_mom(mass, pos, vel, 10*r_half)  # angular momentum vector
    rot_mat = vector_rotation_matrix(L_vec, [0,0,1])  # [0,0,1] for face-on
    new_pos = np.dot(rot_mat, pos.T)
    new_vel = np.dot(rot_mat, vel.T)
    new_pos, new_vel = new_pos.T, new_vel.T
    
    Quen_posG[str(i)]=new_pos
    Quen_velG[str(i)]=new_vel
    
        
    H1,X1,Y1= np.histogram2d(new_pos[:,0], new_pos[:,1], bins=nbins, range=rng)

    H1[H1<=0]=np.nan

    axs[3].imshow(np.log10(H1.T) ,vmin=np.log10(np.nanmin(H1))*1.05,cmap=CMAP,interpolation = 'spline16',interpolation_stage='rgba',origin='lower',extent=[-extent*r_half,extent*r_half,-extent*r_half,extent*r_half])
    axs[3].set_xlabel('x [kpc]')
    axs[3].set_ylabel('y [kpc]')
    
    H1,X1,Y1= np.histogram2d(new_pos[:,0], new_pos[:,2], bins=nbins, range=rng)

    H1[H1<=0]=np.nan

    axs[4].imshow(np.log10(H1.T) ,vmin=np.log10(np.nanmin(H1))*1.05,cmap=CMAP,interpolation = 'spline16',interpolation_stage='rgba',origin='lower',extent=[-extent*r_half,extent*r_half,-extent*r_half,extent*r_half])
    axs[4].set_xlabel('x [kpc]')
    axs[4].set_ylabel('z [kpc]')
    
    H1,X1,Y1= np.histogram2d(new_pos[:,1], new_pos[:,2], bins=nbins, range=rng)

    H1[H1<=0]=np.nan

    IM=axs[5].imshow(np.log10(H1.T),vmin=np.log10(np.nanmin(H1))*1.05,cmap=CMAP,interpolation = 'spline16',interpolation_stage='rgba',origin='lower',extent=[-extent*r_half,extent*r_half,-extent*r_half,extent*r_half])
    axs[5].set_xlabel('y [kpc]') 
    axs[5].set_ylabel('z [kpc]')
    
    plt.colorbar(IM,cax=cax2,ticks=[np.log10(np.nanmin(H1))*1.07,np.log10(np.nanmax(H1))*0.95])
    cax2.set_ylabel('Density', rotation=270)
    cax2.set_yticklabels(['low','high'])
    plt.show()
 
    print('\n')
    
print('\n\n\n')

538464.1418081385
Gas: Count:  494611

To few particles showing instead DM
Gas: Count:  0

To few particles showing instead DM
Gas: Count:  4102

To few particles showing instead DM
Gas: Count:  0

To few particles showing instead DM
Gas: Count:  0

To few particles showing instead DM
Gas: Count:  7

To few particles showing instead DM
Gas: Count:  486

To few particles showing instead DM
Gas: Count:  0

To few particles showing instead DM
Gas: Count:  0

To few particles showing instead DM
Gas: Count:  0

SFR_pos={}
SFR_vel={}

SFR_posG={}
SFR_velG={}

DM_or_gas_SF={}

for i in Save_Galaxies['50']['SF']:
    
    
    fig= plt.figure(figsize=(18,9))
    gs = gridspec.GridSpec( 3,7,                             # Number of axis y,x
                            height_ratios = [1,0.23,1],             # relatives ratios of heigh
                            width_ratios  = [1,0.1,1,0.1,1,0.1,0.08], # relatives ratio of with
                            left  = 0.1,      # Space to the edge of left   from the nearest axis
                            right = 0.95,     # Space to the edge of right  from the nearest axis
                            bottom= 0.15,     # Space to the edge of bottom from the nearest axis
                            top   = 0.97,     # Space to the edge of top    from the nearest axis
                            wspace= 0.1,      # Space horizontal between each of the axis
                            hspace= 0.0)      # Space vertical   between each of the axis

    ax1  = fig.add_subplot(gs[0,0])
    ax2  = fig.add_subplot(gs[0,2])
    ax3  = fig.add_subplot(gs[0,4])
    cax1 = fig.add_subplot(gs[0,6])
    
    ax4  = fig.add_subplot(gs[2,0])
    ax5  = fig.add_subplot(gs[2,2])
    ax6  = fig.add_subplot(gs[2,4])
    cax2 = fig.add_subplot(gs[2,6])
    

    axs=[ax1,ax2,ax3,ax4,ax5,ax6]
    
    
    fig.suptitle(f'Galaxy {i} / Star Forming',y=1.05)
    #################################################################################
    
    subhalo_id, r_half, com, cov = get_basic_info_of_subhalo(ID=i)
    fields = ['Coordinates', 'Velocities', 'Masses']

    st = il.snapshot.loadSubhalo(basePath,snapNum,subhalo_id,'stars',fields=fields)

    mass = st['Masses']
    pos  = st['Coordinates']
    vel  = st['Velocities']

    pos = pos - com  # shifting co-ordinate system to center of subhalo
    vel = vel - cov


    # plotting three projection of subhalo after rotation
    #fist rotate subhalo 
    L_vec = cal_specific_ang_mom(mass, pos, vel, 10*r_half)  # angular momentum vector
    rot_mat = vector_rotation_matrix(L_vec, [0,0,1])  # [0,0,1] for face-on
    new_pos = np.dot(rot_mat, pos.T)
    new_vel = np.dot(rot_mat, vel.T)
    new_pos, new_vel = new_pos.T, new_vel.T
    
    SFR_pos[str(i)]=new_pos
    SFR_vel[str(i)]=new_vel
    
        
    extent = 5  # in unit of half mass radius
    rng = [[-extent*r_half,extent*r_half],[-extent*r_half,extent*r_half]]
    nbins = int((rng[0][1] - rng[0][0])/0.12)
    
    H1,X1,Y1= np.histogram2d(new_pos[:,0], new_pos[:,1], bins=nbins, range=rng)

    H1[H1<=0]=np.nan

    axs[0].imshow(np.log10(H1.T), vmin=np.log10(3) ,cmap='gist_ncar',interpolation = 'spline36',interpolation_stage='rgba',origin='lower',extent=[-extent*r_half,extent*r_half,-extent*r_half,extent*r_half])
    axs[0].set_xlabel('x [kpc]')
    axs[0].set_ylabel('y [kpc]')
    
    
    H1,X1,Y1= np.histogram2d(new_pos[:,0], new_pos[:,2], bins=nbins, range=rng)

    H1[H1<=0]=np.nan

    axs[1].imshow(np.log10(H1.T), vmin=np.log10(3) ,cmap='gist_ncar',interpolation = 'spline36',interpolation_stage='rgba',origin='lower',extent=[-extent*r_half,extent*r_half,-extent*r_half,extent*r_half])
    axs[1].set_xlabel('x [kpc]')
    axs[1].set_ylabel('z [kpc]')
    
    H1,X1,Y1= np.histogram2d(new_pos[:,1], new_pos[:,2], bins=nbins, range=rng)

    H1[H1<=0]=np.nan

    IM=axs[2].imshow(np.log10(H1.T), vmin=np.log10(3) ,cmap='gist_ncar',interpolation = 'spline36',interpolation_stage='rgba',origin='lower',extent=[-extent*r_half,extent*r_half,-extent*r_half,extent*r_half])
    axs[2].set_xlabel('y [kpc]')
    axs[2].set_ylabel('z [kpc]')
    
    plt.colorbar(IM,cax=cax1,ticks=[np.log10(3)*1.05,np.log10(np.nanmax(H1))*0.95])
    cax1.set_ylabel('Density', rotation=270)
    cax1.set_title('Stars')
    cax1.set_yticklabels(['low','high'])
    
    ####################################################################################
    ####################################################################################
    ##############    GAS or DM   ######################################################
    ####################################################################################
    ####################################################################################
    
    try:
        st = il.snapshot.loadSubhalo(basePath,snapNum,subhalo_id,'gas',fields=fields)
        
        if st['count']<100:
            raise Exception('To few particles showing instead DM') 
        
        mass = st['Masses']
        pos  = st['Coordinates']
        vel  = st['Velocities']
        cax2.set_title('Gas')
        
        print('Gas: Count: ',st['count'])
        
        CMAP = 'inferno'
        DM_or_gas_SF[str(i)]='Gas'
        
    except Exception as e:
        print(e)
        
        print('Gas: Count: ',st['count'])
        
        fields = ['Coordinates', 'Velocities']
        st = il.snapshot.loadSubhalo(basePath,snapNum,subhalo_id,'dm',fields=fields)
    
        pos  = st['Coordinates']
        vel  = st['Velocities']
        mass = np.repeat(MASS_DM,len(pos))
        cax2.set_title('DM')
        CMAP = 'viridis'
        DM_or_gas_SF[str(i)]='Dark Matter'

    pos = pos - com  # shifting co-ordinate system to center of subhalo
    vel = vel - cov


    # plotting three projection of subhalo after rotation
    #fist rotate subhalo 
    L_vec = cal_specific_ang_mom(mass, pos, vel, 10*r_half)  # angular momentum vector
    rot_mat = vector_rotation_matrix(L_vec, [0,0,1])  # [0,0,1] for face-on
    new_pos = np.dot(rot_mat, pos.T)
    new_vel = np.dot(rot_mat, vel.T)
    new_pos, new_vel = new_pos.T, new_vel.T
    
    SFR_posG[str(i)]=new_pos
    SFR_velG[str(i)]=new_vel
    
        
    H1,X1,Y1= np.histogram2d(new_pos[:,0], new_pos[:,1], bins=nbins, range=rng)

    H1[H1<=0]=np.nan

    axs[3].imshow(np.log10(H1.T) ,vmin=np.log10(np.nanmin(H1))*1.05,cmap=CMAP,interpolation = 'spline16',interpolation_stage='rgba',origin='lower',extent=[-extent*r_half,extent*r_half,-extent*r_half,extent*r_half])
    axs[3].set_xlabel('x [kpc]')
    axs[3].set_ylabel('y [kpc]')
    
    H1,X1,Y1= np.histogram2d(new_pos[:,0], new_pos[:,2], bins=nbins, range=rng)

    H1[H1<=0]=np.nan

    axs[4].imshow(np.log10(H1.T) ,vmin=np.log10(np.nanmin(H1))*1.05,cmap=CMAP,interpolation = 'spline16',interpolation_stage='rgba',origin='lower',extent=[-extent*r_half,extent*r_half,-extent*r_half,extent*r_half])
    axs[4].set_xlabel('x [kpc]')
    axs[4].set_ylabel('z [kpc]')
    
    H1,X1,Y1= np.histogram2d(new_pos[:,1], new_pos[:,2], bins=nbins, range=rng)

    H1[H1<=0]=np.nan

    IM=axs[5].imshow(np.log10(H1.T),vmin=np.log10(np.nanmin(H1))*1.05,cmap=CMAP,interpolation = 'spline16',interpolation_stage='rgba',origin='lower',extent=[-extent*r_half,extent*r_half,-extent*r_half,extent*r_half])
    axs[5].set_xlabel('y [kpc]') 
    axs[5].set_ylabel('z [kpc]')
    
    plt.colorbar(IM,cax=cax2,ticks=[np.log10(np.nanmin(H1))*1.07,np.log10(np.nanmax(H1))*0.95])
    cax2.set_ylabel('Density', rotation=270)
    cax2.set_yticklabels(['low','high'])
    plt.show()
 
    print('\n')

print('\n\n\n')

Gas: Count:  351436

Gas: Count:  71575

Gas: Count:  349190

Gas: Count:  110659

Gas: Count:  358181

Gas: Count:  67241

Gas: Count:  15699

Gas: Count:  65494

Gas: Count:  5464

Gas: Count:  126201

2. Compute their stellar and gas density profiles between 0 and their stellar half-mass radius.

## To save the data of the galaxies
# np.save('SFR_pos.npy', SFR_pos)  
# np.save('SFR_vel.npy', SFR_vel)  
# np.save('Quen_pos.npy', Quen_pos)  
# np.save('Quen_vel.npy', Quen_vel)  

N=100

fig= plt.figure(figsize=(18,9))
gs = gridspec.GridSpec( 3,9,                             # Number of axis y,x
                        height_ratios = [1,0.2,1],             # relatives ratios of heigh
                        width_ratios  = [1,0.05,1,0.05,1,0.05,1,0.05,1], # relatives ratio of with
                        left  = 0.1,      # Space to the edge of left   from the nearest axis
                        right = 0.95,     # Space to the edge of right  from the nearest axis
                        bottom= 0.15,     # Space to the edge of bottom from the nearest axis
                        top   = 0.97,     # Space to the edge of top    from the nearest axis
                        wspace= 0.1,      # Space horizontal between each of the axis
                        hspace= 0.0)      # Space vertical   between each of the axis

ax1  = fig.add_subplot(gs[0,0])
ax2  = fig.add_subplot(gs[0,2],sharey=ax1)
ax3  = fig.add_subplot(gs[0,4],sharey=ax1)
ax4  = fig.add_subplot(gs[0,6],sharey=ax1)
ax5  = fig.add_subplot(gs[0,8],sharey=ax1)

ax6  = fig.add_subplot(gs[2,0],sharey=ax1)
ax7  = fig.add_subplot(gs[2,2],sharey=ax1)
ax8  = fig.add_subplot(gs[2,4],sharey=ax1)
ax9  = fig.add_subplot(gs[2,6],sharey=ax1)
ax10  = fig.add_subplot(gs[2,8],sharey=ax1)

AX=[ax1,ax2,ax3,ax4,ax5,ax6,ax7,ax8,ax9,ax10]

AX[0].set_yscale('log')

COLOR={'Stars':'blue','Dark Matter':'springgreen','Gas':'red'}

for j,i in enumerate(Quen_pos.keys()):
    subhalo_id, r_half, com, cov = get_basic_info_of_subhalo(ID=int(i))
    
    #####################################################
    # Stars
    #####################################################
    
    Dato  = Quen_pos[i]
    
    R_part= (Dato[:,0]**2+Dato[:,1]**2+Dato[:,2]**2)**0.5
    R_part= R_part[(R_part<=r_half)]
    
    H,E = np.histogram(R_part,bins=N)
    
    W= ((4/3)*np.pi*(E[1:]**3)) - ((4/3)*np.pi*(E[:-1]**3))
    #V= (4/3)*np.pi*r_half**3
    
    AX[j].stairs(H/W,E/r_half,label='Stars',color=COLOR['Stars']) # Divided by the volume enclosed , normilized to the half radius
    AX[j].set_title(f'Galaxy {i}/ Quenched',fontsize=10)
    
    #####################################################
    # Gas or DM
    #####################################################
    
    Dato  = Quen_posG[i]
    
    TYPE = DM_or_gas_Q[str(i)]
    
    R_part= (Dato[:,0]**2+Dato[:,1]**2+Dato[:,2]**2)**0.5
    R_part= R_part[(R_part<=r_half)]
    
    H,E = np.histogram(R_part,bins=N)
    
    W= ((4/3)*np.pi*(E[1:]**3)) - ((4/3)*np.pi*(E[:-1]**3))
    #V= (4/3)*np.pi*r_half**3
    
    AX[j].stairs(H/W,E/r_half,label=TYPE,color=COLOR[TYPE]) # Divided by the volume enclosed , normilized to the half radius
    LG=AX[j].legend(loc='upper right',fontsize=9)
    LG.set_title(f'$R_{{half}}: {r_half:.3f}$',prop={'size':10})
    
    

    
# Erasing extra y labels
plt.setp(ax2.get_yticklabels(), visible=False)
plt.setp(ax3.get_yticklabels(), visible=False)
plt.setp(ax4.get_yticklabels(), visible=False)
plt.setp(ax5.get_yticklabels(), visible=False)
plt.setp(ax7.get_yticklabels(), visible=False)
plt.setp(ax8.get_yticklabels(), visible=False)
plt.setp(ax9.get_yticklabels(), visible=False)
plt.setp(ax10.get_yticklabels(), visible=False)

plt.setp(ax1.get_xticklabels(), visible=False)
plt.setp(ax2.get_xticklabels(), visible=False)
plt.setp(ax3.get_xticklabels(), visible=False)
plt.setp(ax4.get_xticklabels(), visible=False)
plt.setp(ax5.get_xticklabels(), visible=False)

fig.supylabel('Log ([$N_{particules}/kpc^3$])',fontsize=18)
fig.supxlabel('Log ($R/R_{Half}$)',fontsize=18)

plt.show()
print('\n')  

fig= plt.figure(figsize=(18,9))
gs = gridspec.GridSpec( 3,9,                             # Number of axis y,x
                        height_ratios = [1,0.2,1],             # relatives ratios of heigh
                        width_ratios  = [1,0.05,1,0.05,1,0.05,1,0.05,1], # relatives ratio of with
                        left  = 0.1,      # Space to the edge of left   from the nearest axis
                        right = 0.95,     # Space to the edge of right  from the nearest axis
                        bottom= 0.15,     # Space to the edge of bottom from the nearest axis
                        top   = 0.97,     # Space to the edge of top    from the nearest axis
                        wspace= 0.1,      # Space horizontal between each of the axis
                        hspace= 0.0)      # Space vertical   between each of the axis

ax1  = fig.add_subplot(gs[0,0])
ax2  = fig.add_subplot(gs[0,2],sharey=ax1)
ax3  = fig.add_subplot(gs[0,4],sharey=ax1)
ax4  = fig.add_subplot(gs[0,6],sharey=ax1)
ax5  = fig.add_subplot(gs[0,8],sharey=ax1)

ax6  = fig.add_subplot(gs[2,0],sharey=ax1)
ax7  = fig.add_subplot(gs[2,2],sharey=ax1)
ax8  = fig.add_subplot(gs[2,4],sharey=ax1)
ax9  = fig.add_subplot(gs[2,6],sharey=ax1)
ax10  = fig.add_subplot(gs[2,8],sharey=ax1)

AX=[ax1,ax2,ax3,ax4,ax5,ax6,ax7,ax8,ax9,ax10]

AX[0].set_yscale('log')

for j,i in enumerate(SFR_pos.keys()):
    subhalo_id, r_half, com, cov = get_basic_info_of_subhalo(ID=int(i))
    
    #####################################################
    # Stars
    #####################################################
    
    Dato  =  SFR_pos[i]
    
    R_part= (Dato[:,0]**2+Dato[:,1]**2+Dato[:,2]**2)**0.5
    R_part= R_part[(R_part<=r_half)]
    
    H,E = np.histogram(R_part,bins=N)
    
    W= ((4/3)*np.pi*(E[1:]**3)) - ((4/3)*np.pi*(E[:-1]**3))
    #V= (4/3)*np.pi*r_half**3
    
    AX[j].stairs(H/W,E/r_half,label='Stars',color=COLOR['Stars'])
    AX[j].set_title(f'Galaxy {i}/ Star Forming',fontsize=10)
    
    #####################################################
    # Gas or DM
    #####################################################
    
    Dato  = SFR_posG[i]
    
    TYPE = DM_or_gas_SF[str(i)]
    
    R_part= (Dato[:,0]**2+Dato[:,1]**2+Dato[:,2]**2)**0.5
    R_part= R_part[(R_part<=r_half)]
    
    H,E = np.histogram(R_part,bins=N)
    
    W= ((4/3)*np.pi*(E[1:]**3)) - ((4/3)*np.pi*(E[:-1]**3))
    #V= (4/3)*np.pi*r_half**3
    
    AX[j].stairs(H/W,E/r_half,label=TYPE,color=COLOR[TYPE]) # Divided by the volume enclosed , normilized to the half radius
    LG=AX[j].legend(loc='upper right',fontsize=9)
    LG.set_title(f'$R_{{half}}: {r_half:.3f}$',prop={'size':10})


    
# Erasing extra y labels
plt.setp(ax2.get_yticklabels(), visible=False)
plt.setp(ax3.get_yticklabels(), visible=False)
plt.setp(ax4.get_yticklabels(), visible=False)
plt.setp(ax5.get_yticklabels(), visible=False)
plt.setp(ax7.get_yticklabels(), visible=False)
plt.setp(ax8.get_yticklabels(), visible=False)
plt.setp(ax9.get_yticklabels(), visible=False)
plt.setp(ax10.get_yticklabels(), visible=False)

plt.setp(ax1.get_xticklabels(), visible=False)
plt.setp(ax2.get_xticklabels(), visible=False)
plt.setp(ax3.get_xticklabels(), visible=False)
plt.setp(ax4.get_xticklabels(), visible=False)
plt.setp(ax5.get_xticklabels(), visible=False)

fig.supylabel('Log ([$N_{particules}/kpc^3$])',fontsize=18)
fig.supxlabel('Log ($R/R_{Half}$)',fontsize=18)

plt.show()
print('\n')  
    

3. Select the dark matter particles of the 2 most massive dark matter haloes in TNG100. Compute their mass density profiles between the softening length and the virial radius.

# import astropy.units as u

# p_crit=((1.9e-29)*(h**2))*u.g*(u.cm**-3)
# p_crit=p_crit.to(u.Msun*(u.kpc**-3))
# p_crit

# p_crit=p_crit.value

We can obtain the softening lenght from the webpage of TNG

(https://www.tng-project.org/about/)

i        = 100
basePath = '/home/tnguser/sims.TNG/TNG'+str(i)+'-1/output/'

FIELDS=['GroupFirstSub','GroupMass','GroupMassType','Group_M_Crit200','Group_R_Crit200']
HHH=il.groupcat.loadHalos(basePath,99,fields=FIELDS)


GFSub = HHH['GroupFirstSub']

GMass = HHH['GroupMass']

Mtype = HHH['GroupMassType']

M200  = HHH['Group_M_Crit200']

R200  = HHH['Group_R_Crit200']


# https://www.tng-project.org/data/downloads/TNG100-1/
h       = 0.6774
MASS_DM = (0.000505574296436975e10)/h # M sun 

Halos = {}
ddd   = GMass.copy()

if PRINT:
    print('# Most Massive v/s Group First Sub')

for i in range(3):
    a=ddd.argmax()
    if PRINT:
        print(f' {a}\t{GFSub[i]}')
    ddd[a]=0
    #pos      = il.groupcat.loadSingle(basePath,99,subhaloID=GroupFirstSub[i])
    subhalo_id, r_half, com, cov = get_basic_info_of_subhalo(Sim='100',ID=a)
    
    pos      = il.snapshot.loadSubhalo(basePath,99,a,'dm',fields=['Coordinates'])
    pos      = pos-com
    R_part   = ((pos[:,0]**2+pos[:,1]**2+pos[:,2]**2)**0.5)/h
    Halos[a] = R_part
    
del ddd


if PRINT:
    print(Halos.keys())

# Most Massive v/s Group First Sub
 0	0
 1	17185
 2	31342
dict_keys([0, 1, 2])

The most massive Halos of TNG 100-1

fig= plt.figure(figsize=(18,9))
gs = gridspec.GridSpec( 1,5,                             # Number of axis y,x
                        height_ratios = [1],             # relatives ratios of heigh
                        width_ratios  = [1,0.05,1,0.05,1], # relatives ratio of with
                        left  = 0.1,      # Space to the edge of left   from the nearest axis
                        right = 0.95,     # Space to the edge of right  from the nearest axis
                        bottom= 0.15,     # Space to the edge of bottom from the nearest axis
                        top   = 0.97,     # Space to the edge of top    from the nearest axis
                        wspace= 0.1,      # Space horizontal between each of the axis
                        hspace= 0.0)      # Space vertical   between each of the axis

ax1  = fig.add_subplot(gs[0,0])
ax2  = fig.add_subplot(gs[0,2])#,sharey=ax1)
ax3  = fig.add_subplot(gs[0,4])#,sharey=ax1)


AX=[ax1,ax2,ax3]


N=100

XY_data={}

for j,i in enumerate(Halos.keys()):
    
    R_part = Halos[i]                                  # Get the data stored
    R_trim = R_part[((0.74<R_part)&(R_part<R200[i]))]  # Cut the data in the range asked
    R_trim = np.log10(R_trim)                          # in log space
    
    H,E = np.histogram(R_trim,bins=N)
    
    E = 10**E
    W = ((4/3)*np.pi*(E[1:]**3)) - ((4/3)*np.pi*(E[:-1]**3)) # Volume in each step 
    Y = (H/W)*MASS_DM
    
    XY_data[i]=[mid(E),Y]
    
    AX[j].stairs(Y,E/R200[i]
                ,label=f'Halo {i}')
    
    #AX[j].axhline(200*p_crit,color='red')
    AX[j].set_yscale('log')
    AX[j].set_xscale('log')
    AX[j].legend(title=f'$R_{{200}}$= {R200[i]:.2f}',loc='upper right')

    
fig.supylabel('Log ([M$_{\odot}/kpc^3$])',fontsize=18)
fig.supxlabel('Log (R/$R_{200}$)',fontsize=18)
plt.show()
print('\n\n\n')

4. Try to see if you can fit each measured radial density function with a Navarro-Frenk-White (NFW) density profile.

$$\rho(r,\rho_0,R_s)=\frac{\rho_0}{\frac{r}{R_s}(1+\frac{r}{R_s})^2}$$

# We will use the scipy package in order to get the best fit
from scipy.optimize import curve_fit

def NFW(x,p,r): # Navarro-Frenk-White (NFW) profile.
    X=(x/r)
    return np.log10(p/(X*((1+X)**2))) # We do the fit in Log space
    

fig= plt.figure(figsize=(18,9))
gs = gridspec.GridSpec( 1,5,                             # Number of axis y,x
                        height_ratios = [1],             # relatives ratios of heigh
                        width_ratios  = [1,0.05,1,0.05,1], # relatives ratio of with
                        left  = 0.1,      # Space to the edge of left   from the nearest axis
                        right = 0.95,     # Space to the edge of right  from the nearest axis
                        bottom= 0.15,     # Space to the edge of bottom from the nearest axis
                        top   = 0.97,     # Space to the edge of top    from the nearest axis
                        wspace= 0.1,      # Space horizontal between each of the axis
                        hspace= 0.0)      # Space vertical   between each of the axis

ax1  = fig.add_subplot(gs[0,0])
ax2  = fig.add_subplot(gs[0,2])#,sharey=ax1)
ax3  = fig.add_subplot(gs[0,4])#,sharey=ax1)


AX=[ax1,ax2,ax3]

AAA=[[1e6,300],[1e6,150],[1e7,40]] #Init guess (something reasonable)

for j,i in enumerate(XY_data.keys()):

    xdata=(XY_data[i][0])
    ydata=np.log10(XY_data[i][1]) # We do the fit in Log space

    A=AAA[j]
    
    popt, pcov = curve_fit(NFW, xdata, ydata,maxfev=10000,p0=A,bounds=(0.74, [1e16, 1e3]))
    
    if PRINT:
        print(f'For Halo {i} we have a $rho_0$= {popt[0]:.3f},\t and a $R_s$={popt[1]:.3f}')
        
    AX[j].scatter(xdata/R200[i]  ,ydata                     -np.log10(popt[0]) ,label='Data', marker='o',facecolor='none',edgecolor='r')
    AX[j].plot(   xdata/R200[i]  ,NFW(xdata,popt[0],popt[1])-np.log10(popt[0]) ,label='Best Fit')
    AX[j].plot(   xdata/R200[i]  ,NFW(xdata,A[0],A[1])      -np.log10(popt[0]) ,label='My eye Fit')
    # plt.yscale('log')
    AX[j].set_xscale('log')
    AX[j].set_title(f'Halo {i}')
    LG=AX[j].legend(loc='upper right',fontsize=9)
    LG.set_title('From Best Fit\n'+rf'Log ($\rho_0$)= {np.log10(popt[0]):.3f}'+f'\n$R_s$         = {popt[1]:.2f}',prop={'size':10})
    
fig.supylabel(r'Log ($\rho/\rho_0$)',fontsize=18)
fig.supxlabel(r'Log ($R/R_0$)',fontsize=18)
print()

For Halo 0 we have a $rho_0$= 552491.892,	 and a $R_s$=365.993
For Halo 1 we have a $rho_0$= 4189014.063,	 and a $R_s$=69.198
For Halo 2 we have a $rho_0$= 48541611.218,	 and a $R_s$=14.726

print('\n\n\nBy: Ian Baeza')

By: Ian Baeza